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ABSTRACT 

Thermal instability (TI) can strongly affect the structure and dynamics of the interstellar medium (ISM) in 
the Milky Way and other disk galaxies. Thermal conduction plays an important role in the TI by stabilizing 
small scales and limiting the size of the smallest condensates. In the magnetized ISM, however, heat is con- 
ducted anisotropic ally (primarily along magnetic field lines). We investigate the effects of anisotropic thermal 
conduction on the nonlinear regime of the TI by performing two-dimensional magnetohydrodynamic simula- 
tions. We present models with magnetic fields of different initial geometries and strengths, and compare them 
to hydrodynamic models with isotropic conduction. We find anisotropic conduction does not significantly alter 
the overall density and temperature statistics in the saturated state of the TI. However, it can strongly affect the 
shapes and sizes of cold clouds formed by the TI. For example, for uniform initial fields long filaments of cold 
gas are produced that are reminiscent of some observed HI clouds. For initially tangled fields, such filaments 
are not produced. We also show that anisotropic conduction suppresses turbulence generated by evaporative 
flows from the surfaces of cold blobs, which may have implications for mechanisms for driving turbulence in 
the ISM. 

Subject headings: galaxies: ISM - instabiUties - ISM: kinematics and dynamics 



1. INTRODUCTION 

The thermal instability (TI) plays an important role in con- 
trolling several different aspects of the interstellar medium 
(ISM) and star formation. For example, it has been inv oked to 
explain the observed multiphase str ucture of the ISM dParkei 
[1953; SpitziilUll [Field 1965; Fi eld. Goldsmith. & Habing 
[1969). The linear stage of the TI in astrophysical gases 
was first studied in detail by Field (1965). He identified 
three unstable modes: an isobaric and two isentropic modes. 
In the nonlinear regime, the isobaric mode produces con- 
densations that ar e fundamental to the classical tw o-phase 
model of the ISM (iField. Goldsmifli. & H abing' 1969^. as well 
as the extension to a three-pha se model by Cox & SmithI 
(11974 : iMckee & Ostrikei ^ (1911). The TI also regulates the 
mass flow between the different components of the ISM, 
and therefore affects the star formation rate (Chieze 1987). 
Thus, it is important to investigate the role of TI in deter- 
mining the distribution of density, temperature, and other 
physical variables in the multiphase ISM. For this reason, 
a variety of authors have studied the Hnear and nonlinear 
1 stages of the TI using numeri cal hydrodynamic si mulations 
jHennebelle & Perault) 119991 lKr itsuk & Norman ! I2002a)|^ 
2004; Vazquez-Semadeni, Gazol, & Scalo 2000; Gazol et al. 
2001; Sanchez-Salcedo, Vazquez-Semadeni & Gazol 2002; 
Piontek & Ostriker 2004, 2005; Kovama & Inutsuka 2004; 
Vazquez-Semadeni et al. 20071: iKim. Kim. & Ostrikeri 120081: 
Inoue & Inutsuka 2008). 

Thermal conduction is important to include in studies of 
the TI for two reasons. First, it suppresses the growth rate 
at small scales, in fact isobaric perturbations with wavelength 
smaller than the Field length Xf ( Fieldlll965) do not grow at 
all. Second, it produces evaporation from the surfaces of cold 
dense fragments, and the interaction of the evaporative flows 
can induce turbulence. 

Including thermal conduction is essential for numerical 
studies of the TI. Without explicit thermal conduction, per- 
turbations at the grid scale grow fastest, and may even- 



tuafly come to dominate the dynamics. For this reason, 
iKovama & Inutsukal (120041) concluded that numerical stud- 
ies of the TI must satisfy a "Field criterion", that is the 
Field length must be resolved by at least a few cells to pre- 
vent artificial fragmentation at the grid scale, and to avoid 
the results being dominated by grid noise. Satisfying the 
Field cr iterion requires including ex plicit thermal conduction 
dPionte k & Ostriker 2004, 2005; K ovama & Inutsukal |2004 
'Brandenburg, Korpi & Mee 2007). To highlight one exam- 
ple, Piontek & Ostriker (2004, 2005) studied the interaction 
of the TI and the magnetorotational instability (MRI) in disks 
including isotropic thermal conduction; they found the MRI 
could drive turbulence and fragmentation in the diffuse ISM 
at amplitudes consistent with HI observations. 

To date most studies of the TI with conduction have as- 
sumed the conductivity is isotropic. However, in a magne- 
tized plasma, electrons can flow more freely along magnetic 
field lines t han across the m, leading to anisotropic transport 
coefficients (ISpitzerlll962h . The degree of anisotropy is mea- 
sured by the ratio of the electron gyro radius to the mean free 
path between collisions. For the warm medium, where typi- 
cally T = 1500 K, n = 2 cm^, and B = 1/iG, the Columb mean 
free path is Amfp ^ 10'° cm, while the electron gyro radius 
is rg = 10^ cm. Thus, in this medium the thermal conduction 
should be highly anisotropic, and primarily along magnetic 
field lines. 

The implications of anisotropic transport terms on the dy- 
namics o f astrophysical plasmas has begun to be explored 
recently ('Balbus"2001t ISharma. Hammett & Ouataerti 120031: 
iParrish & Stone 200^. For example, in stratified atmo- 
spheres, anisotropic conduction can result in the magnetother- 
mal and he at-flux bu oyancy in stabilities in the int racluster 
medium dParrish & Stonej 120081: ISharma et al.ll2009l) . In ad- 
dition, it can have effects on th e evolution of supernova 
remnants (iBalsara. TiUev & Howkll2 008) and. on magnetized 
spherical accretion flows Tsii arma. Ouataert & Stonal2008h . 
Recently, .Sharma. Parrish & Ouataert (2010.) studied the TI 
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in the intracluster medium, including heating by cosmic rays. 
Their primary result is that with anisotropic thermal conduc- 
tion, the TI could produce filaments of cold gas orientated 
along magnetic field lines. 

In this study, we perform two-dimensional numerical hy- 
drodynamic and magnetohydrodynamic simulations of TI in 
the ISM that include the anisotropic heat conduction. The pur- 
pose of this paper is to investigate the effects of anisotropic 
thermal conduction on the structure and dynamics of inter- 
stellar medium, and to compare the results to models which 
include only isotropic conduction. Since the geometry and 
strength of magnetic field in different regions of the ISM may 
vary, we perform simulations of the TI with various mag- 
netic field strengths and two different initial geometries for 
the field. 

Ideally, the numerical studies of the thermal instability 
should include both isotropic and anisotropic conduction 
with the temperature dependency since the thermal instabil- 
ity develops density distribution with distinct peaks in the 
cold/dense phase dominated by isotropic conduction and in 
the diffuse/hot phase with higher ionization fraction and thus 
more affected by anisotropic conduction. However, includ- 
ing both isotropic and anisotropic conduction with a realis- 
tic temperature dependency is challenging, since it would de- 
crease the Field length, and require higher resolution. In this 
study, we restrict our exploration to the simulations only with 
isotropic conduction or ones only with anisotropic conduction 
to better understand the effects of anisotropic conduction. 

This paper is organized as follows: our numerical methods 
and code tests are summarized in § |2] In § [51 we first dis- 
cuss the effects of the conduction rate and resolution on the 
TI in hydrodynamical simulations, and then use these results 
to choose a specific set of model parameters for our simula- 
tions. Results from calculations of the nonlinear regime of 
the TI with anisotropic conduction are presented in § |4] We 
summarize and discuss our results in § [S] 

2. NUMERICAL MODEL 

We solve the equations of ideal MHD with additional terms 
for radiative cooling, heating, and heat conduction 

|^ + V-[pv] = 0, (1) 



dpy 



+ V- 



PVV-BB+ \p+%- ] I-cr 



dE „ 



E + P + — ) v-B(B-v) + Q-v-cr 



= 0, (2) 
= -pA(3) 



dB 
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where the symbols have their usual meaning. The total energy 
density E is given as 



E = 



1 



+ -nv +■ , 

7-1 2' 2' 



(5) 



where 7 is the ratio of specific heats (assumed to be 5/3 
throughout this paper). Here, a is the viscous stress tensor 



cr,7 = V 



dxj dxi ) 3 



(6) 



where 77 is the coefficient of viscosity and summation over 
repeated indices is implied. In equation[3] Q denotes the heat 



flux and C = pK(p,T)-T is the net cooling function, where 
A(p. T) and V are the cooling and heating rates per unit mass 
respectively. These equations are written in units such that the 
magnetic permeability pB= 1 ■ 
The heating and cooling rates are 



pmur = 2x 10"^^ ergs 



(7) 
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where we have adopte d the functional fit to the IS M cooling 
rates as suggested by iKoyama & Inutsukal (l2002lPl For the 
net cooling curve adopted in this work, the transition temper- 
atures that define the warm (phase "F", T > Tmax), interme- 
diate (phase "G", r^in < T < T^ax), and cold (phase "H", at 
T < T^in) phases are T^ax = 5012 K and r^in = 185 K respec- 
tively. 

In the case of isotropic conduction, the heat flux Q is given 
by 

Q = Qiso = -/Cvr, (9) 

while in the case of anisotropic conduction, the heat flux is 
given by 

Q = Qani,so=-/Cbb-Vr. (10) 

where b is a unit vector in the direction of the magnetic field, 
and K, is the conductivity. For t he temperature ranges cons id- 
ered here, the latter is given by (iParkeril 1 953t ISpitzerilT962h 



/C = 2.5 X 10^^ r'/^erg cm"' K"' s" 



(11) 



The relation between the dynamic viscosity coefficient 77 
and the thermal conductivity K, can be characterized by the 
Prandtl number: 

1 kB V 



Pr = 



7 - 1 niu K, 



(12) 



The Prandtl number for a neutral monoatomic gas is Pr=2/3. 
For most of the simulations presented in this paper, we include 
explicit viscosity with a constant coefficient of viscosity that 
corresponds to a Prandtl number Pr = 2/3. 

All of the numerical results presented in this paper were 
calculated using Athena. Details of the algorithms imple- 
mented in Athena are documented in Gardiner & Stone ( 20051 
I2OO8I) and ex tensive tests of the MHD algorithms in Athena 
are shown in IStone et akl (l2008h . In this work, we use the 
Roe approximate Riemann solver and the directionally unsplit 
CTU (corner transport upwind) integration scheme. Optically 
thin cooling was added to the integrators without using oper- 
ator splitting. To circumvent timestep constraints in regions 
of strong cooling, we limit the minimum temperature to the 
equilibrium value where heating balances cooling. This al- 
lows large time steps while avoiding the large errors or over- 
shoots encountered with backward Euler or Crank-Nicholson 
implicit differencing. 

We solve the equations in a rectangular domain with peri- 
odic boundary conditions and no gravity. Our domain spans a 

* Equation [s] contains the corrections pointed out by 
INagashima. Inutsuka & Kovamal j200fl) 
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Fig. 1. — Numerically measured growth rates of the TI for different wave 
numbers and with a constant conductivity /C = 10^ erg cm"' K"' s"' plotted 
as points, overlaid with the theoretical growth rates for conductivities of /C = 
[0.0, 10^, 10*, 10', 10^] erg cm"' K"' s"'. The excellent agreement between 
the numerical and analytic growth rates verifies our numerical algorithms. 

region of size 2L x L, where L is varied between different cal- 
culations. We adopt an initial pressure of P/k = 3000 K cm"^ 
and density of n = 2.0 cm"^. The equilibrium pressure corre- 
sponding to this density is P/fc= 2896.6 K cm"^, so our initial 
condition is in a slowly cooling state. The mean mass per 
particle in units of the mass of atomic hydrogen is /i = 1 .27, 
representing 10% He abundance by number To initiate the 
thermal instability (TI), we add random pressure perturbations 
with a maximum amplitude of 0.1%. For our assumed initial 
temperature, the fiducial value of the conductivity is 

/Co = 9.68 X lOVg cm"' K"' s"' (13) 

which follows from equation [TTI 

In order to test the implementation of the heating, cooling 
and conduction terms in our code, we have performed one- 
dimensional simulations of the TI, and compared the numer- 
ically measured growth rates with the theoretical prediction 
iFieldl (11965 ). The initial conditions for this test are a medium 
at rest with constant density and pressure of n = 2.0 cm""* and 
P/k = 3000 K cm"""*, and isotropic conduction with coefficient 
K, = lO^erg s"'cm"'K"'. The box size was L = 100 pc and 
the grid contained 1024 zones. We initialize the models with 
eigen modes of the instability by imposing sinusoidal fluctu- 
ations of amplitude 1 percent and with wave number k, and 
measure the growth rate of the density for different values 
of k. Figure [1] compares the numerical growth rates from 
our simulations along with the theoretical values using four 
values of the thermal conductivity /C. Clearly there is good 
agreement between the analytic and numerical growth rates; 
the numerical method reproduces the theoretical value of the 
growth rate to better than 4 percent in all cases. 

For any non-zero value of the conductivity, the thermal dif- 
fusivity becomes comparable to the heating and cooling term 
at a critical wave length referre d to as the Field length, given 
approximately bv lFieldl (Il965h 

A. = 2.y||. (14) 

In the case of no conduction /C = 0, the growth rate of the 
TI is largest at the smallest scales (A = 0). For this reason. 



TABLE 1 

Models with isotropic conduction. 



Conductivity 




Box size (pc) 




A.r(pc) 




0.0 




4.0x2.0 




0.00195 




Ka 


0.12/0.006 


4.0 X 2.0 


33.3 


0.00195 


61.5/3 


4/Co 


0.24/0.012 


8.0 X 4.0 


33.3 


0.0039 


61.5/3 


16/Co 


0.48 / 0.024 


16.0x8.0 


33.3 


0.0078 


61.5/3 




0.96 / 0.048 


32.0 X 16.0 


33.3 


0.0156 


61.5/3 



"Field length in parsec as given by Eauation ll4l 

''The number of Field lengths along the longest axis of the simulation box 
(2L). 

■^Field number, i.e., the number of zones in a Field length. 

'Koyama & Inutsuka (2004) have shown that the inclusion of 
explicit thermal conduction is necessary in studies of the TI in 
order to damp growth at the grid sca le and ensure the fastest 
growing modes of the TI lFieldl (11965 ) are resolved by at least 
three cells. Without explicit conduction, grid noise will be 
amplified until it dominates the solution. We further inves- 
tigate the consequences of simulating the TI without explicit 
thermal conduction in the next section. 

3. TI WITH ISOTROPIC THERMAL CONDUCTION 

In the following subsections, we investigate the role of 
varying the amplitude of the conductivity, size of the com- 
putational domain, and grid resolution using models with 
isotropic conduction. This allows us to decide on an opti- 
mal set of model parameters before considering the effect of 
anisotropic conduction in the next section. 

3.1. Effect of varying the conductivity 

With our adopted initial condition of density and temper- 
ature and with the standard level of thermal conduction JCq, 
the Field length \f is about 0.12 pc. Resolving this length 
in a computational domain that spans hundreds of parsecs is 
very demanding. Thus, in most previous work, the thermal 
conductivity has either been increased (b y about a factor of 
100) in order to resolve the Field length (iPiqntek & Ostrikerl 
lIOOllIOOl 'Branden burg. Korpi & Meel2"007ir or the compu- 
tational domain has been chosen to be very small (for example 
0.3-4.8pc) lKovama & Inutsuka (2004, 2006) so that the stan- 
dard value /Co can be used. In this section, we study whether 
the properties of the nonlinear regime of the TI depend on the 
value of the conductivity adopted. 

We present the results of five inviscid hydrodynamical mod- 
els which include isotropic conduction with values for the 
conductivity of /C = [0.0,/Co,4/Co, 16/Co, 64/Co]- Table 1 sum- 
marizes the properties of the simulations. 

Each model begins from a uniform homogeneous medium 
at rest with the fiducial density and pressure. Each uses the 
same numerical resolution, 2048 x 1024 cells, but the size of 
the domain L is varied so that the Field length is resolved, even 
in the nonlinear regime. Since the Field length is proportional 
to \/T and inversely proportional to density n, as density in- 
homogeneities are amplified by the TI the local value of the 
Field length varies, i.e., the local Field length in the regions of 
hot gas increases and that in regions of cold gas decreases. By 
the end of the simulation, we find the Field length in the cold 
gas is ^ 20 times smaller than in the initial conditions. To 
fully resolve the Field length throughout the simulation, grid 
cells must be smaller than the Field length in the cold phase. 
We calculate the minimum value of the Field length Af min for 
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Mass-weighted density PDFs 
19.90Myr 



Fig. 2. — Density snapshots at f 200 Myr for the inviscid TI simulations 
without thermal conduction (top) and with isotropic uniform conduction of 
/C = /Co (bottom). A logarithmic color scale from n = 0.28 cm"^ to n= 170 
cm"' is used. 

cold gas and keep the Field number ric = Xp/Sx, greater than 
3 throughout our simulations. 

We choose the physical size of the domain L to be large 
enough so that it spans a large number of Field le ngths, sim- 
ilar to the value used in previous w ork ( Piontek & Ostriken 
l200ll200irBrandenburg. Korpi & Meai2007 ). Defining ha^ 
as the number of Field lengths along the longest axis of the 
simulation box (2L), we vary L so that for each value of the 
conductivity, n\p is fixed. We have also performed a large 
number of simulations with the same conductivity K, but dif- 
ferent values of L to confirm that the statistical properties of 
the density and the level of kinetic energy in the nonlinear 
stage of the TI are not affected by changing L (provided that 
n\p is large). Thus, we are confident we can compare runs 
with different physical sizes to isolate the effect of varying 
the conductivity on the TI. 

Figure |2] plots the density at f ~ 200 Myr in two models: 
a calculation without conduction (K, = 0, top) and one with 
K, = /Co (bottom). Density structures formed by the TI are sig- 
nificantly different between the two models. In the case of 
no thermal conduction (top panel), very small dense and cold 
fragments begin to form at the beginning of the simulation. 
These fragments are unresolved (1-2 cells in size). As time 
proceeds they begin to merge, and by the end of the calcula- 
tion at 200 Myr, they dominate the results. On the other hand, 
for the simulation with /C = /Co (bottom panel), resolved den- 
sity structures begin to emerge at about 20 Myr, and continue 
to grow into a network of filaments that forms by about 25 
Myr By ^ 33 Myr, dense cold clouds in thermal equilibrium 
form. The shape and size of the cold clouds evolves due to 
evaporation and winds driven from their surfaces by thermal 
conduction. The smallest clouds evident in the density figure 
have a short lifetime (less then ^ 2 Myr), as they are quickly 
dissolved by evaporation, and continuously reformed by the 
evolving density structures. 

Figure [3] shows the mass-weighted density PDFs for the 




-1.0-0,5 0.0 0.5 1.0 1.5 2.0 -1.0-0.5 0.0 0.5 1.0 1.5 2.0 2.5 

loq (n) 

Fig. 3. — Mass weighted density PDFs for inviscid TI simulations with 
different values of an isotropic conductivity axt ^ 20, 25, 33, 200 Myr. 
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Fig. 4. — Mass weighted mean velocity (top) and velocity dispersion (bot- 
tom) of inviscid TI models with different values of an isotropic conductivity 
as a function of time. 

five models with different level of thermal conductivity. At 
t = 20 Myr, all the gas in the models with thermal conduction 
are in the unstable regime, however in the model without con- 
duction the density is already segregated into two phases. Not 
until f = 33 Myr do the models with conduction form den- 
sity distributions that are clearly separated into two distinct 
phases. After this time, the shape and size of dense clouds 
evolve due to condensation, evaporation, and thermal winds 
as shown in Figure|2] so that there are always some fraction of 
the gas in the interface regions that is in the thermally unstable 
regime. The mass of gas in the unstable regime increases as 
the value of the conductivity increases. Moreover, the highest 
density also increases with the value of the conductivity. 

In Figure |4] we plot the mass-weighted mean velocity, de- 
fined as V = ^ < pv^ > / < >, and the velocity dispersion 
for each of the five models. In each case the velocity dis- 
persion increases rapidly by Tl-induced turbulence. Without 
thermal conduction (orange solid line), the velocity disper- 
sion and average velocity attain their maximum values by f 
20Myr. In contrast, in the models with thermal conduction, 
the velocity dispersion and the average velocity keep increas- 
ing because of the continuous supply of unstable gas produced 
by thermal evaporation from the surfaces of den se clouds as 
discussed in .Inoue. Inutsuka & Kovamai (120061) . Thus, the 
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Fig. 5. — The saturated velocity as a function of isotropic tliermal conduc- 
tivity. 

model without thermal conduction (orange solid line) has the 
smallest values of both the average velocity and velocity dis- 
persion, while the model that includes the highest level of con- 
ductivity (red dotted line) has the largest. 

The mass-weighted mean velocity calculated at f > 50 Myr 
as a function of the thermal conductivity is shown in Figure 
|5] Clearly the amplitude of the turbulence induced by thermal 
instability increases in proportion to the value of the thermal 
conductivity, at least for /C > /Co- 

We draw two conclusions from this study. First, 

as has alread y been sho wn in many recent s tudies 

jKovama & Inutsuka" l2004t iPiontek & Ostriked 12004), we 
conclude it is critically important to include explicit thermal 
conduction in simulations of the TI. Figure |2] dramatically 
demonstrates how calculations without conduction lead to un- 
resolved fragments that do not evaporate or drive turbulence. 
Figure [3] shows the substantial differences between the PDF 
of the density in models that include conduction, in compar- 
ison to a model that does not. Any study of Tl-driven tur- 
bulence in the ISM that does not include conduction will be 
dominated by errors seeded by the grid. Second, we conclude 
it is important to adopt the appropriate value for the conduc- 
tivity. Adopting a larger value in order to resolve the Field 
length results in overestimating the level of turbulence driven 
by evaporative flows, and modifies the PDF of the density. 

3.2. Effect of numerical resolution and viscosity 

Based on the results from the previous subsection, in the re- 
mainder of this paper we present models in which the conduc- 
tivity /C = /Co- The corresponding Field l ength is Af = 0.12pc. 
Based on one dimensional simulations, iKovama & Inutsuka 
(|2004) have shown that the Field number = ^F.min/Sx 
should be greater than 3 to achieve convergence in the number 
of clouds formed by the TI, and the maximum Mach number 
at late time in their simulations. It is important to confirm 
that convergence is achieved if the Field length is resolved 
in our two dimensional study. Therefore we have run invis- 
cid simulations of the TI with isotropic conduction with three 
different numerical resolutions: 1024 x 512 (low resolution), 
2048 X 1024 (standai-d resolution), and 4096 x 2048 (high res- 
olution). All models are computed in a domain of size 2L x L, 
where L = 2pc. We have also run models with a constant coef- 
ficient of dynamic viscosity corresponding to a Prandtl num- 
ber Pr = 2/3 at all three resolutions. 

Figure |6] shows the density PDFs at f ^ 150 Myr both for 
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Fig. 6. — Mass weiglited density PDFs at / ~ 150Myr in models computed 
witli different resolutions. (Top) Density PDFs of inviscid models. (Bottom) 
Density PDFs of models with explicit viscosity and Pr = 2/3. 
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Fig. 7. — Mass weighted velocity dispersion as a function of time at dif- 
ferent numerical resolutions for inviscid calculations (top), and calculations 
with explicit viscosity and Pr = 2/3. 

inviscid models and models with viscosity at different reso- 
lutions. The PDF of the inviscid calculations (top panel) are 
resolution dependent, especially for the maximum density and 
the amount of mass in the unstable regime. The highest reso- 
lution inviscid model shows the largest value of the maximum 
density reached in cold gas, and has a much larger amount of 
matter in the unstable regime compared to lower resolutions. 
The latter is a reflection of the fact that matter in the unstable 
regime is contained in dynamical regimes that are driven by 
thermal evaporation, and the properties of these evaporative 
flows are affected by viscosity. The bottom panel of Figure |6] 
shows that both the maximum value of the density and the 
amount of matter in the unstable regime converge with reso- 
lution if explicit viscosity is included. In this case, the mass- 
weighted fraction of matter in the unstable regime is 5.8, 8.2, 
9.0 % for the low, standard, and high resolution models re- 
spectively. 
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TABLE 2 

Models with anisotropic conduction 



iVHJQCl 


l^r\n Hi 1 r*l"i rMi 

\^ Ul 1 LI LH^ LI tJ 1 1 


V laHJally 


CTnf*ti H 1 H 


P 


1 


hydro/Isotropic 


inviscid 






2 


mhd/Anisotropic 


inviscid 


uniform 




3 


mhd/Anisotropic 


inviscid 


uniform 


1 


4 


mhd/Anisotropic 


inviscid 


tangled 




5 


mhd/Anisotropic 


inviscid 


tangled 


1 


6 


hydro/Isotropic 


Pr = 2/3 






7 


mhd/Anisotropic 


Pr = 2/3 


uniform 


lO*^ 


8 


mhd/Anisotropic 


Pr = 2/3 


uniform 


1 


9 


mhd/Anisotropic 


Pr = 2/3 


tangled 


10* 


10 


mhd/Anisotropic 


Pr = 2/3 


tangled 


1 



Figure |7] shows the mass- weighted velocity dispersion as a 
function of time at different numerical resolutions for mod- 
els both with and without explicit viscosity . In all the mod- 
els, the velocity dispersion starts to rapidly increase at around 
f ^ 20 Myr, when dense structures begin to form (see also 
Figures|2]and[3j. While the velocity dispersion of the inviscid 
models continues to increase until f ^ 50 Myr, in the mod- 
els with explicit viscosity it remains relatively constant after 
turbulence initially develops. In either case, the velocity dis- 
persion saturates after f ^ 50 Myr, however the mean velocity 
in the saturated state of the inviscid models is almost 5 times 
larger than those that include explicit viscosity. Specifically, 
for f > 50 Myr the mass-weighted mean velocity in the in- 
viscid models with low, standard, and high resolution is 0. 10, 
0. 19, 0.24 km/s respectively, while models with viscosity have 
mean velocities that are only 0.017, 0.037, 0.043 km/s. This 
difference is likely due to the fact that the inviscid models 
have 2-3 times more matter in the unstable regime, which 
drives correspondingly more powerful evaporative flows. The 
highest resolution model with viscosity has only a 15% larger 
mean velocity in the saturated state than that of the standard 
resolution model, which is further evidence all quantities are 
converging. 

Based on these convergence tests, we conclude that some 
properties of the nonlinear regime of the TI with isotropic 
conduction are still resolution dependent even if the Field 
length is well resolved if explicit viscosity is not included. 
For this reason, in studies with anisotropic conduction pre- 
sented in the next section, we will present results both with 
and without viscosity. 

4. TI WITH ANISOTROPIC THERMAL CONDUCTION 

In this section, we investigate the effect of anisotropic con- 
duction on the development of the TI by performing magneto- 
hydrodynamics (MHD) simulations and comparing the results 
to hydrodynamic simulations of the TI with isotropic conduc- 
tion. All of the simulations start with our fiducial density 
and pressure, use a conductivity JC = /Co, and are computed 
in a domain of size 4 x 2 pc with our standard resolution of 
2048 X 1024. This means the Field length is resolved even in 
the cold phase with at least 3 cells. 

4.1. Model parameters 

We have performed a series of simulations in order to ex- 
plore both the effect of varying the magnetic field geometry 
and strength on the nonlinear outcome of the TI. To explore 
the effect of geometry, we have computed models with an ini- 
tially uniform magnetic field inclined along the diagonal of 




Fig. 8. — Density snapshots at f ~ 200 Myr in TI simulations with isotropic 
conduction (top), anisotropic conduction with weak uniform field (second 
panel), strong uniform field (third panel), and strong tangled field (bottom) 
respectively. The arrow shows the initial direction of the magnetic field. 

the domain (to eliminate potential artifacts that might be pro- 
duced by aligning the field with the grid), as well as a ran- 
dom tangled field generated using the technique described 
below. In both cases, we study two different initial mag- 
netic field strengths corresponding to /3 = Pgas/Pmag = [lO*", 1]. 
For our fiducial initial conditions, these field strengths are 
B = [0.00323, 3.23] AfG. 

To generate the tangled magnetic field, we initialize a vector 
potential A = (A^ . A,,, A,) which has only one non-zero compo- 
nent that is given by a Fourier power spectrum of the form 

A. = |(5At| cxr"/^ (15) 
with amphtudes that follow a Gaussian random distribution 
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Fig. 9. — Mass (black solid line) and volume (blue dotted line) weighted 
density PDFs for models 6. 7, 8 and 10 at f ~ 24 and 200 Myr. The percent- 
ages of gas in each phase by mass are noted in black, and the proportions by 
volume are shown in blue. 

for all wavenumbers k in the range 2tt/{L/2) < k < 2tt/6x. 
The magnetic field is then B = V x A, and is normalized so 
that the volume averaged magnetic energy gives the appropri- 
ate value of /3. 

Table 2 summarizes the properties of our runs. They in- 
clude two models with isotropic conduction (one with and one 
without viscosity), as well as 8 runs with anisotropic conduc- 
tion with different magnetic field strengths and geometry. We 
study the TI with anisotropic conduction using both models 
that are inviscid, and models that include explicit viscosity 
with Pr = 2/3. 

In this paper, we present four of our simulations: a 
hydrodynamic model with the isotropic conduction in an 
unmagnetized medium (Model 6) and MHD models with 
the anisotropic conduction in magnetized medium (Model 
7, 8, and 10). For detailed study on the formation of 
clouds by the TI with isotropic conduction or without ther- 
mal conduction in a magnetized medium, we refer to re- 
cent pa pers: Hennebelle et al. (2008); Heitsch et al (200^; 
llnoue & Inutsukal (120091) : iBanerjee et al.1 (120091) . 

4.2. Evolution of the density 

Figure [8] shows snapshots of the density at f ^ 200 Myr 
from four of our simulations: a hydrodynamic model with 
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Fig. 10. — Mass (black solid line) and volume (blue dotted Hne) weighted 
temperature PDFs for models 6, 7, 8 and 10 at f 24 and 200 Myr The 
percentages of gas in each phase by mass are shown. 

isotropic conduction, MHD models with anisotropic conduc- 
tion and both a weak and strong uniform field, and an MHD 
model with anisotropic conduction and a strong tangled field. 
In each case, the TI results in cold regions that are more than 
an order of magnitude denser than their surroundings, how- 
ever, the distribution and evolution of the cloud blobs are sig- 
nificantly different in each case. 

For the hydrodynamic model with isotropic conduction 
(top), a network of filaments forms at about 25 Myr connect- 
ing the regions of highest density. By ^ 33 Myr, dense cold 
clouds begin to form. Clouds with a radius smaller than the 
Field length Af are destroyed by thermal conduction. Evap- 
orative flows form around larger cold clouds, they are evi- 
dent as the "fuzzy" edges of dense regions in Figure [8] The 
evaporative flows cause the dense clouds and filaments both 
to merge and fragment in random motions. 

For the MHD model with anisotropic conduction and ei- 
ther a weak or strong uniform magnetic field, the regions of 
cooled gas tend to form thin filaments parallel to the direction 
of the field lines. In the linear regime, the preference for fil- 
amentation along the fields can be explained by the variation 
of the Field length with angle with respect to the direction of 
field. Perpendicular to the field, the Field length is very small, 
so very short wavelength perturbations are unstable. Parallel 
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to the field, the Field length is much longer, so only longer 
wavelength modes grow. As the perturbations grow nonlin- 
ear, conduction tends to enforce isothermality along magnetic 
field lines, leading to long filaments. In the case of strong 
fields, only motions along field lines are allowed, and more- 
over magnetic pressure provides some support in dense re- 
gions, resulting in fragments aligned with the field lines as 
evident in the third panel of Figure [8] 

For the model with a strong tangled magnetic field (bottom 
panel), spherical dense clouds quickly develop at the centers 
of regions of closed field lines (the local minimum and maxi- 
mum of the vector potential). The topology of the closed field 
lines in these regions means they become thermally isolated 
from the rest of the domain, and therefore neither evaporate 
nor evolve further Small clouds originally formed in regions 
of open field lines eventually evaporate along the magnetic 
field lines and disappear. 

4.3. Density and temperature statistics 

Figure |9] shows the mass and volume weighted density 
PDFs at both 24 (left column) and 200 (right column) Myr for 
the four different models shown in Figure [HJ from top to bot- 
tom the panels are for isotropic conduction, anisotropic con- 
duction with weak and strong uniform field, and anisotropic 
conduction with a strong tangled field respectively. As was 
shown in section 13.11 thermal conduction can reduce the 
growth rate of the TI. However, models with anisotropic con- 
duction develop density distributions with distinct peaks in the 
dense (cold) and diffuse (hot) phases earlier than the model 
with isotropic conduction. Similarly, by the end of the sim- 
ulation, the hydrodynamic model with isotropic conduction 
has a larger amounts of mass in the unstable regime com- 
pared to the models with anisotropic conduction. This in- 
terpret this result as due to the suppression of evaporative 
flows from the surfaces of dense clouds by the magnetic field 
and anisotropic conduction. As shown earlier, material in the 
unstable regime is contained in evaporative flows (the fuzzy 
edges of the clouds seen in the top panel of Figure [8]). With 
anisotropic conduction, thermal conduction can only occur 
from the surfaces of the clouds where the direction of the 
magnetic field is normal. For the elongated and filamentary 
clouds formed with anisotropic conduction, only a very small 
fraction of its surface area is subject to evaporation, result- 
ing in very little unstable gas. For the same reason, the frac- 
tion of high density gas (cold phase) is slightly larger in the 
anisotropic models compared to that of the isotropic case. 

Probably the most striking feature of the density PDFs are 
the distinct peaks representing the two phase medium at 200 
Myr. The mass and volume weighted temperature PDFs of 
these four models shown in Figure [TOl reflect this structure: 
most of the matter is either in the cold or hot phases. Very lit- 
tle is contained in the unstable regime at intermediate temper- 
atures, and the fraction of gas in the unstable regime decreases 
in the anisotropic conduction models. 

4.4. Evolution of the kinetic energy 

As shown in section lTTl the TI drives evaporative flows that 
produce turbulence. Figure [TT] compares the mass weighted 
mean velocity and velocity dispersion in models that include 
isotropic conduction, anisotropic conduction with a weak uni- 
form field or strong uniform field, and a strong tangled field. 
The evolution of these quantities is a direct measure of the 
strength of evaporative flows. 
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Fig. 11. — Mass weighted mean velocity {top) and velocity dispersion 
(bottom) as a function of time from models with hydrodynamic conduction 
(/3 = inf), 1, 6, 7, 8. 
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Fig. 12. — Contour plot of the magnetic field strength (color) and magnetic 
field lines at ; ~ 200 Myr in TI simulation with anisotropic conduction and 
with weak uniform field (Model 7). 

Except in the case of the tangled field, the mean velocity 
and dispersion grow initially, until dense clouds form and the 
TI saturates around 25 Myr. In the case of the strong tangled 
field, there are unbalanced Lorentz forces in the initial condi- 
tions that cause motions immediately. These decay away as 
the TI saturates. Thus, the time evolution of the velocity and 
its dispersion is not a good measure of TI driven turbulence in 
this case. For the other two models with anisotropic conduc- 
tion (the strong and weak uniform field cases), the saturation 
amplitude of the mean velocity and dispersion at late times 
are about 5-10 times lower in the models with anisotropic 
conduction. This is clear evidence that evaporative flows from 
the surfaces of the dense clouds is suppressed by the magnetic 
field and anisotropic conduction. This is because evaporation 
can only occur when the field line is normal to the interface. 
For most of the surface of long filaments, the field is parallel 
to the interface. This, only a very small region at the end of 
the filament can produce evaporative flows. 

4.5. Amplification of the magnetic field 

In Figure [12] we show an image of the magnetic field 
strength overlaid with line segments to show the direction 
of the magnetic field lines at f ^ 200 Myr in a TI simula- 
tion with anisotropic conduction and with weak uniform field 
(Model 7). The magnetic field is initially uniform with a 
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strength of 0.00323 fiG. As material accumulates into the fil- 
aments by TI, the magnetic fields are compressed, resulting 
in strong fields that are parallel to the direction of the fila- 
ment. Outside of the filaments, the magnetic field is folded 
and twisted by vortical fluid motions. We find the compres- 
sion and twisting can amplify t he field up to 0.3 /i G, i.e., by a 
factor of 100, as shown in Ino ue. Inutsuka & Kova ma (2007). 
In full three-dimensions, the geometrical compression during 
the formation of filaments is likely to be larger than the two- 
dimensional case studied here, thus the factor of 100 ampli- 
fication may be a lower bound. Note the regions where the 
magnetic field is strongest trace the cold hydrogen filaments. 

5. SUMMARY 

In this work, we have studied the nonlinear regime of the 
thermal instability (TI) and the effect of anisotropic thermal 
conduction by performing two-dimensional hydrodynamical 
and magnetohydrodynamic simulations incorporating radia- 
tive cooling and heating. Our main results can be summarized 
as the following: 

1 . As found in previous studies, it is crucial to include ex- 
plicit thermal conduction in numerical studies of the TI so 
that the Field length is resolved, in order to prevent artificial 
fragmentation driven by numerical noise at the grid scale (e.g. 
Figure |2]). 

2. The amplitude of the thermal conductivity controls the 
rate of evaporation from the surfaces of dense clouds formed 
by the TI, and therefore strongly affects the amplitude of tur- 
bulent motions induced by the TI (e.g. Figure|5]l. 

3. Even when the Field length is resolved, explicit vis- 
cosity must be included to obtain numerical convergence of 
some quantities, for example the amplitude of turbulent mo- 
tions driven by evaporative flows (e.g. Figure|7]). 

4. Although the statistics of the density and temperature are 
not strongly affected by anisotropic conduction, the geome- 
try of structures formed by the TI are quite different. With 
anisotropic conduction and a uniform magnetic field, the TI 
saturates as long thin filaments of dense gas aligned with the 
field. In a tangled field, spherical clouds are formed in regions 
of closed field lines (e.g. Figure|8]l. 

5. The combination of anisotropic conduction and MHD 
strongly suppresses the rate of evaporation of cold gas from 
the surfaces of dense structure in regions where the field is 
parallel to the interface. This reduces the amplitude of turbu- 
lence driven by the TI (e.g. Figure [TTTl. 

These results have a number of implications for observa- 
tions of cold neutral gas in the ISM. In particular, the thin 
filaments along the magnetic field in the weak uniform mag- 
netic field case agree well with rec ent observations of the 
Riegel -Crutcher cloud conducted by iMcClure-Griffiths et alj 
( 120061) . This neutral hydrogen (HI) cloud lies on the edge 
of the Local Bubble (Crutcher & Lien 1984) filled with a hot 
and diffuse gas where the anisotropic conduction can be ex- 
pected. McClure-Griffiths et al. (2006 ) found a network of 
dozens of hairlike filaments of cold hydrogen with widths of 
less than ~ 0.1 pc and up to 17 pc long. They also have found 
that the filaments are aligned with the magnetic field of the 
cloud which agrees well with our results. They also calcu- 
lated the magnetic field strength by using the Chandrasekhar- 
Fermi method, finding ^ 60/iG. We find that compression and 
twisting of field during the formation of filaments by the TI 
with anisotropic conduction can amplify it by up to a factor 
of 100 in two-dimensions. In three-dimensions, the amplifi- 
cation is Ukely to be larger. This may be enough to explain 



the observed field; 

Recently, Sha rma. Parrish & Ouataerll (1201 Oh have reported 
a study of the TI with anisotropic thermal conduction in the 
hot X-ray emitting plasma in clusters of galaxies. The heat- 
ing and cooling processes in this regime are very different 
than those in the ISM studied here (equations|7]and[8]respec- 
tively), in particular there are no stable phases in the cooling 
curve they adopt, so that magnetic pressure sets the only limit 
on how cold and dense the gas becomes. Nonetheless, the 
structure of density condensations in this case are very simi- 
lar to our results: filaments of cold gas along magnetic field 
lines. 

Our results show that in the case of anisotropic conduction, 
the geometry of the magnetic field with respect to the interface 
between cold and hot phases is very important. Only if the 
field is normal to the interface can thermal conduction drive 
evaporation and outflows. The inclusion of anisotropic con- 
duction could have important implications for the structure 
of interf aces, and the interpretation of observations of these 
regions ("In oue. Inutsuka & Kovamal 120061 : IStone & Zweibell 
Il009, 2010). 

There are a number of limitations to our work that should 
be addressed in future in vestigations. Firstly, the ISM is 
highly turbulent dHeiles & Troiandl l20()3h . with typical tur- 
bulen t velocities approximately 7 kms~' dHeiles & Trolandl 
120031: iMohan. Dwarakanath & SrinivasanI 120041) . i.e. more 
than 100 times larger than those produced by the TI with 
isotropic conduction and viscosity. In the traditional pic- 
ture of the ISM, sup ernova driven turbul e nce leads to a hot, 
diffus e third phase (Cox & SmithI 119 741: iMckee & Ostrikeil 
Il977h . More recently, (Piontek & Ostri keii 120041) have co"i> 
sidered the interaction of turbulence driven by the MRI and 
TI. We have not considered the effect of externally forced tur- 
bulence on the TI in this work, however this would be a pro- 
ductive direction for study in the future. 

Secondly, in this work we have adopted a constant con- 
ductivity JC. However, in rea l ity the conductivity is a func- 
tion of the temperature iParke J (Il953h . so that the rate of ther- 
mal conduction decreases as the gas cools. The amplitude of 
the conductivity can vary by 2 orders of magnitude between 
warm and cold phases at the end of our typical simulation 
at t ^ 100 Myr. As discussed in Section [TT] the value of 
the conductivity can affect the rate of evaporation from dense 
clouds, and therefore the kinetic energy of turbulence driven 
by the TI, thus assuming a constant value may alter the result. 
Using a realistic temperature dependent conductivity is chal- 
lenging, since it would decrease the minimum value of Field 
length and require much higher resolution. Since the geome- 
try of the magnetic field limits the amount of evaporation to 
a very small surface area at the ends of the filaments in the 
case of anisotropic conduction, we do not expect the use of 
a more realistic conductivity to produce qualitative changes 
in our result. Nonetheless, more realistic studies which use 
temperature dependent conductivities would be fruitful. 

Finally, this study has considered flows in only two dimen- 
sions. In full three dimensions, the amplification of the mag- 
netic field due to geometrical compression into filaments may 
be larger, and the turbulence driven by evaporative flows may 
be of a different character Fully three-dimensional simula- 
tions of the TI with anisotropic conduction would also be in- 
teresting for future studies. 
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